NPSOR-92-012 


NAVAL 


POSTGRADUATE  SCHOOL 

Monterey,  California 


DTIC 

ELECTS 

JUL 15 1992 


MODELING  AND  STATISTICAL  ANALYSIS  OF 
MEDAKA  BIOASSAY  DATA 


Donald  P.  Gaver  and  Patricia  A.  Jacobs 


May  1992 


Approved  for  public  release;  distribution  is  unlimited. 
Prepared  for: 

Army  Biomedical  Research  and  Development  Laboratory 
Ft.  Detrick  MD 


NAVAL  POSTGRADUATE  SCHOOL, 
MONTEREY,  CALIFORNIA 


Rear  Admiral  R.  W.  West,  Jr.  Harrison  Shull 

Superintendent  Provost 


This  report  was  funded  by  the  Army  Biomedical  Research  and  Development 
Laboratory,  Ft  Detrick,  MD. 

This  report  was  prepared  by: 


DONALD  P.  GAVER 


Professor, 


Department  of  Operations  Research 


Reviewed  by: 


PETER  PURDUE 
Professor  and  Chairman 
Department  of  Operations  Research 


PATRICIA  A.  JACOfiS 
Professor, 

Department  of  Operations  Research 
Released  by: 


Paul  j.  marto 

Dean  of  Research 


|gWJilllimW8SlHL«fcl»LWKa«ilMaMq 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMBNo  0704-0183 


m :  i  a  aws  i  oci  Mil :  1 i  u«j  i  a  if  a  i  (• : 


UNCLASSIFIED 


•]  i  w  ?  f/iV/j  i  ViUii  I  m«]ana:Mz«fl 


Approved  for  public  release;  distribution  is 
unlimited. 


NPSOR-92-012 


I  m  i  a  y  ^ :  1 :  i  ,T,  i  k'  [ei«] :  ray  k’  i  i  m  f  ■  ra  *•:  a  a  w  ai  ■ 


Naval  Postgraduate  School 


ity,  State,  and 


Monterey,  CA  93943 

8a  NAME  OF  FUNDlNd/S  PON  SORING 
ORGANIZATION 

Army  Biomedical 
Research  and 
Development  Laborator 


8c.  AD  DRESS  ( City,  State,  and  ZIP  Code ) 

Ft.  Detrick,  MD 


I 


(If  applicable) 

OR 


8b  OFFICE  SYMBOL 
(II  applicable) 


iJPJAiMaAMfiiMiiifieMirayjiwiiiraj 


Naval  Postgraduate  School 


ny.  State,  and  AH 


jEgflagssai 


Monterey,  CA  93943-5006 


PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 


10.  SOURCE  OF  FUNDING  NUMBERS 


PROGRAM  PROJECT  TASK 

ELEMENT  NO.  NO.  NO 


WORK  UNIT 
ACCESSION  NO 


unty  Classification ) 

Modeling  and  Statistical  Analysis  of  Medaka  Bioassay  Data 


liv 


Donald  P.  Gaver  and  Patricia  A.  Jacobs 


Technical  Report 


Ilia— 11  'in  mi  ——i 

GROUP 

SUB-GROUP  1 

w—m—m 

l---  ' 

1 

II:»si»]:»iaH«l3IIA«ir»f 


I  LHMb  ( Continue  on  reverse  it  necessary  and  idenvty  by 

Binomial  Distribution;  censored  data;  generalized  linear 
model;  bootstra 


A  histopathologic  examination  of  tissues  from  Oryzias  latipes  (Japanese  medaka  fish)  was  performed 
to  evaluate  the  carcinogenic  potential  of  tricholoroethylene  (TCE)  in  groundwater.  The  data  were 
reported  by  Experimental  Pathology  Laboratories,  Inc.,  in  a  report  dated  Jan.  19,  1990,  submitted  to  the 
Army  Biomedical  Research  and  Development  Laboratory,  Ft.  Detrick,  MD. 

This  paper  provides  a  brief  statistical  analysis  of  some  aspects  of  those  data.  The  analysis  does  not 
reveal  a  strong  positive  relationship  between  TCE  concentration  over  the  range  considered  and 
probability  (risk  or  hazard)  of  incurring  at  least  one  end  point  manifestation  (here  cystic  degeneration  or 
liver  neoplasm)  in  a  fish.  Uncertainties  in  the  point  estimates  are  assessed  by  bootstrapping.  Both  non- 
parametric  (weak  statistical  assumptions)  and  parametric  (stronger  statistical  assumptions)  analyses  give 
similar  inconclusive  dose-response  indications. 

A  brief  discussion  is  included  of  a  biologically-based  mathematical  model  that  is  likely  to  form  an 
appropriate  basis  for  more  sophisticated  data  analysis. 

One  contribution  of  this  paper  is  to  discuss  and  illustrate  techniques  for  quantitative  analysis  of 
other  similar  data.  The  methods  can  also  be  used  to  assist  in  choosing  un  experimental  design. 


kiHi]  let  I  ziuiviiw  ?  f/AYiWJ  z\ 1 1 1  if •]  |:MH  I 


[7J  UNCIASSIFIEDOJNLIMITED  Q  SAME  AS  RPT.  []  DTIC  USERS  UNCLASSIFIED 


D.  Gaver 


orm  14 


MaaM*ttAis»]i 


Previous  editions  are  obsolete. 

S/N  0102-LF -014-6603 


OR/Gv 


UNCLASSIFIED 


I : :  U'H  f.LS-i  I  d  W.WW  J  f.T*  fi  I H  3 


MODELING  AND  STATISTICAL  ANALYSIS  OF  MEDAKA  BIOASSAY  DATA 
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ABSTRACT 

A  histopathologic  examination  of  tissues  from  Oryzias  latipes  (Japanese  medaka 
fish)  was  performed  to  evaluate  the  carcinogenic  potential  of  tricholoroethylene 
(TCE)  in  groundwater.  The  data  were  reported  by  Experimental  Pathology 
Laboratories,  Inc.,  in  a  report  dated  Jan.  19,  1990,  submitted  to  the  Army  Biomedical 
Research  and  Development  Laboratory,  Ft.  Detrick,  MD. 

This  paper  provides  a  brief  statistical  analysis  of  some  aspects  of  those  data.  The 
analysis  does  not  reveal  a  strong  positive  relationship  between  TCE  concentration 
over  the  range  considered  and  probability  (risk  or  hazard)  of  incurring  at  least  one 
end  point  manifestation  (here  cystic  degeneration  or  liver  neoplasm)  in  a  fish. 
Uncertainties  in  the  point  estimates  are  assessed  by  bootstrapping.  Both  non- 
parametric  (weak  statistical  assumptions)  and  parametric  (stronger  statistical 
assumptions)  analyses  give  similar  inconclusive  dose-response  indications. 

A  brief  discussion  is  included  of  a  biologically-based  mathematical  model  that  is 
likely  to  form  an  appropriate  basis  for  more  sophisticated  data  analysis. 

One  contribution  of  this  paper  is  to  discuss  and  illustrate  techniques  for 
quantitative  analysis  of  other  similar  data.  The  methods  can  also  be  used  to  assist  in 
choosing  an  experimental  design. 


INTRODUCTION 


The  Japanese  medaka  fish  ( Oryzias  latipes )  has  come  to  be  of  great  interest  as  an 
indicator  of  groundwater  toxicity;  see  Van  Beneden  et  al.  (1990),  Gardner,  et  al.  (1990). 

The  Research  Model  Branch,  Health  Effects  Research  Division  of  the  Army 
Biomedical  Research  and  Development  Laboratory,  Ft.  Detrick,  MD,  has  conducted 
extensive  experimentation  with  medaka  so  as  to  test  its  response  to  various  known 
or  suspected  toxic  agents  or  carcinogens.  This  paper  provides  a  statistical  analysis  of 

data  from  such  an  experiment.  Analysis  provides  a  quantitative  and  focussed _ 

perspective  on  the  message  of  the  data  that  usefully  supplements  the  more  usual  I,or 
simple  qualitative  observations. 


Design  of  Experiment 

The  experiment  whose  data  is  analyzed  was  planned  and  conducted  as  follows. 
Eight  (8)  groups  of  medaka  were  treated  as  shown  in  Table  1.  Those  groups  treated 
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with  DEN  received  pretreatment  with  10  mg/1,  of  diethylnitrosamine  for  48  hours  at 
17  days  after  hatch.  The  groups  that  received  TCE  received  various  concentrations  of 
trichloroethylene  (100%,  50%,  25%  and  0%)  on  a  biologically-motivated  scale:  100% 
refers  to  undiluted  groundwater  containing  TCE,  and  50%  and  25%  refer  to 
corresponding  with  pure  water. 


TABLE  1 

TREATMENT  COMBINATIONS  AND  CD  RESPONSES 


#  fish  with  symptom/#  fish 
killed 

Sacrifice  Time 

Group 

DEN 

no  DEN 

%TCE 

3  months 

6  months 

1 

X 

0 

6/25 

4/15 

3 

X 

25 

4/25 

5/13 

5 

X 

50 

2/25 

4/14 

7 

X 

100 

3/25 

3/14 

2 

X 

0 

11/25 

6/12 

4 

X 

25 

4/25 

8/13 

6 

X 

50 

6/25 

5/12 

8 

X 

100 

7/25 

3/8 

The  individual  fish  were  assigned  to  tanks  of  water,  presumably  maintained  at 
standard  temperature,  also  presumably  in  a  random  manner.  There  do  not  appear  to 
have  been  replicate  tanks.  After  three  months  an  interim  sacrifice  was  made  of  25 
fish  in  each  group;  the  number  of  fish  showing  cystic  degeneration  (CD),  after  three 
months/six  months  appear  to  the  left  of  the  slash  (/)  in  the  table.  Thus  Group  1 
contained  six  fish  out  of  25  with  CD  after  three  months,  and  four  out  of  15  after  six 
months;  in  the  latter  case  15  fish  were  exposed  to  the  original  concentration  for  the 
entire  six  months;  this  is  referred  to  as  the  chronic  group.  Another  group,  the  so- 
called  recovery  group,  was  placed  in  pure  water  for  the  second  three  month  period. 
This  group's  response  is  not  analyzed  in  this  paper.  Table  2  reports  the  incidence  of 
liver  neoplasms  for  the  same  fish  in  groups  2,  4,  6,  8  that  were  pretreated  with  DEN. 


TABLE  2 

INCIDENCE  OF  LIVER  NEOPLASMS  IN  MEDAKA  FOR  GROUPS  PRETREATED 

WITH  DEN 


#  fish  with  symptom/#  fish  killed 
Sacrifice  time 

Group 

%  TCE 

3  Months 

6  Months 

2 

0 

0/25 

1/12 

4 

25 

2/25 

3/13 

6 

50 

0/25 

1/12 

8 

100 

2/25 

4/8 

2 


Model-based  Analysis 

The  experimental  outcomes  are  viewed  from  the  following  perspective.  Each 
individual  fish  subjected  to  a  particular  DEN-TCE  treatment  (e.g.,  DEN  =  0, 
TCE  =  50%)  is  initially  thought  of  as  a  member  of  a  population  of  similar  fish.  For 
various  reasons,  including  that  of  genetic  diversity,  the  individual  fish  will  exhibit 
particular  symptoms,  i.e.,  reach  specified  biological  endpoints  such  as  cystic 
degeneration  or  neoplasm  within  specified  organs,  at  widely  different  times.  In 
addition  some  fish  may  die  before  any  such  symptoms  manifest  themselves. 
Consequently  it  is  reasonable  and  natural  to  think  of  the  occurrence  of  a  particular 
endpoint  as  a  probabilistic  (or  random,  or  stochastic)  phenomenon,  much  as  a 
coinflip  or  dice  throw  outcome  is  thought  of,  or  in  the  same  way  that  actuarial 
scientists  regard  human  life  durations  when  proposing  life  insurance  contracts.  That 
is,  let  T,  the  time  to  occurrence  of  a  particular  endpoint,such  as  cystic  degeneration, 
be  a  random  variable,  a  quantity  whose  value  (for  a  particular  fish)  is  determined  by 
sampling  from  a  population  with  a  fixed  distribution  function  that  in  turn  depends 
upon  the  treatment  of  interest  (DEN  and  TCE),  but  also  upon  water  temperature  and 
presence  of  other  elements  in  the  fish  tank,  and  also  individual  fish  traits.  This 
distribution  function  is 

Population  Fraction  of  Fish  with  Symptom 

Times,  T,  less  than  t(t  =  3  months,  6  months)  - 

The  population  parameter  values  that  identify  the  particular  distribution  are  called 
0  =  (01/  #2,  •••/  #p)-  For  instance,  9\  might  be  the  population  mean,  and  62  the 
population  standard  deviation.  Occasionally  the  specific  distribution  used  to  model 
the  variability  in  a  population  of  times  is  normal,  i.e.  its  density  function 


g(t;9)  =  dG  /  dt  =  e  2V  2 


the  familiar  bell-shaped  curve: 


However,  more  frequently  that  data  variability  is  better  described  as  log-normal : 
logarithm  of  T  =  X  has  a  bell-shaped  density  for  the  raw  T-values  tend  to  straggle  off 
to  the  right,  i.e.  the  density  of  T  is  (possibly)  "positively  skewed."  A  simpler  form 
that  may  be  appropriate  is  the  exponential: 
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G(f;A)  =  1  -  e~*t. 

In  this  model  A  =  l/(Mean  time  to  exhibit  symptoms,  e.g.,  CD).  We  will  later  use  this 
exponential  model  in  an  illustrative  analysis.  Still  another  form  that  may  be 
appropriate  for  describing  the  time  to  the  onset  of  a  cancerous  growth  is  the  Weibull 
which  describes  an  increasing,  or  decraesing,  time  of  exposure  effect,  depending 
upon  data  requirements.  Later  we  shall  describe  a  distribution  that  arises  from 
plausible  biological  assumptions,  particularly  when  cancer  is  considered — the 
Moolgavkar  family  of  clonal  expansion  models;  Moolgavkar  et  al.  (1979).  Use  of  the 
latter  "biologically-based"  family  requires  that  at  least  four  parameter  values  be 
estimated  from  data.  In  the  light  of  the  current  experimental  design  such  models  are 
somewhat  difficult  to  identify.  The  exponential  distribution  will  be  used  in  this 
report  to  illustrate  a  parametric  analysis  (one  using  a  specific  assumed  form  for 

atm. 

Use  of  the  Conceptual  Model 

Since  TCE  is  a  toxic  substance  it  might  be  anticipated  that  (a)  the  mean  or 
average  fraction  of  fish  exposed  to  x%  of  TCE  that  exhibit  symptoms  after  t  =  3 
months  (the  first  three  months)  would  increase  with  x  (the  TCE  concentration); 
likewise  for  six  months  of  (chronic)  exposure;  and  (b)  that  the  mean  fraction  of  the 
fish  that  survive  the  first  three  months  that  exhibit  symptoms  in  the  second  three 
months  (between  three  months  and  six  months  of  exposure)  might  increase  over  the 
mean  fraction  in  the  first  three  months.  The  latter  behavior  would  imply  an 
increasing  hazard  property  attributable  to  dosage  with  TCE.  The  increasing  hazard 
property  is  consistent  with  the  idea  that  prolonged  exposure  (to  TCE  here)  increases 
the  chances  that  particular  endpoints  will  occur  as  time  goes  on.  However,  evidence 
for  such  behavior  from  the  current  data  is  not  strong  in  the  light  of  the  uncertainties 
associated  with  sampling  errors  as  assessed  by  bootstrapping. 

NON-PARAMETRIC  ANALYSIS 

Suppose  the  above  general  sampling  model  prevails.  Then  an  estimate  of  the 
probability  that  a  fish  exhibits  a  particular  symptom,  e.g.,  CD,  within  time  t  (=  3 
months)  is 

Estimate  of  Probability  _  - 
That  T  <  t  (=  3  months)  = 

Number  of  Fish  Sacrificed  at  t(  =  3)  that  Exhibit  Symptom  (e.g.,  CD)  ' 
Number  of  Fish  Exposed  (for  t  =  3) 

This  estimate  is  easily  calculated  for  all  treatments;  however  sometimes  it  is  zero. 
This  suggests  that  the  number  of  fish  exposed  is  too  small  to  be  detecta bly  influenced 
by  the  dosage;  in  general  we  might  expect  some  response.  Since  G(3;0)  is  the  so-called 
hazard  associated  with  the  appearance  of  the  particular  symptom  during  the  first 

three  months  of  exposure,  G(3;8)  is  an  estimate  thereof  on  the  basis  of  only  25 
exposed  fish;  for  a  different  25  fish  treated  equivalently  we  generally  anticipate  a 
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different  numerical  value  of  G.  By  re-sai„ fjling  ("bootstrapping)  it  is  possible  to 
appraise  the  sample  variation  in  the  estimate  G:  sample  from  a  binomial 

distribution  with  G  being  the  probability  of  "success"  =  symptom  occurrence  within  t 
=  3,  and  N(=  25),  the  number  exposed,  to  obtain  a  pseudo  or  bootstrapped  sample 
number  of  fish  exhibiting  the  symptom,  and  from  this,  divided  by  N,  a  possible 

sample  value  of  hazard,  i.e.,  Gi.  Repeat  to  get  G2,  again  to  obtain  63,...,  Gb,  where  B 
is  "large."  The  sampling  has  been  repeated  B  =  500  times.  Then  compute 


Variance  G  =  g 


(fib  —  Gb  ) 


2 


and  the  standard  error  of  the  original  estimate,  G,  is  SE[G]  =  V  variance  G  where  Gb 

is  the  sample  mean  of  the  bootstrap  estimates;  Gg  =  F°r  the  first  three- 

month  data  the  above  standard  error  can  actually  be  calculated  directly  (no  re- 

VG(1  -  G) 

— -  ,  but  this  formula  approach  is  not  so  easy 

for  the  second  three-month  period.  Roughly  speaking,  the  true  value  of  G(3)  lies 

within  G-2SE[G]  and  G  +  2SE[G].  So  an  estimate,  and  an  error  estimate,  for  initial 
three  months  hazard  is  obtained.  See  Table  3  for  quoted  estimates  and  standard 
errors  (in  parentheses). 

To  compare  to  the  second  three  months'  hazard  compute 
Estimate  of  Probability 

That  T<t  =  6  Months,  o\ 

Given  that  T  >  f  =  3  Months  =  -^-4 - —  s  C(6,3;0)) 

1-GO;© 

=  Estimate  of  Second  3-Months'  Hazard. 

Notice  that  the  estimates  G(6;0)and  G(3;0)  must  be  obtained  from  different  sets  of 
counts,  just  as  was  done  earlier,  and  consequently  that  there  is  no  guarantee  that 

G(6 ,0)  is  greater  than  G(3 ;6).  Although  no  case  of  such  reversal  occurs  in  the  present 
data,  a  few  reversals  have  occurred  when  resampling  or  bootstrapping  is  done;  in 
such  cases  the  hazard  value  is  set  equal  to  zero.  Standard  errors  of  the  second  3- 

months'  hazards  are  calculated  by  bootstrapping  G(6,2;0)  by  resampling  for  each 
component,  G(6;0)and  G(3;0),  and  combining  as  in  the  formula  above. 

Tables  3  and  4  summarize  the  results  of  the  point  estimates  and  their  standard 
errors.  Table  3  refers  to  CD,  while  Table  4  addresses  neoplasms.  Figures  1  through  5 
graphically  display  the  actual  hazard  sampling  variations  as  assessed  by 

bootstrapping.  Figures  1-3  present  boxplots  of  the  bootstrap  sample  hazards  Gfc(3;0) 
and  Gj,(6,3;0) . 
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The  following  description  of  the  boxplot  is  taken  from  the  documentation  of 
GRAFSTAT,  a  developmental  product  of  IBM  which  the  Naval  Postgraduate  School 
is  using  under  a  test  agreement  with  IBM.  'The  box  portion  of  the  plot  extends  from 
the  lower  quartile  of  the  sample  to  the  upper  quartile.  (The  lower  quartile  is  the 
point  for  which  one  quarter  of  the  sample  lies  below  and  three  quarters  above.  The 
upper  quartile  is  analogous.)  The  line  across  the  center  of  the  box  marks  the  median. 
The  circle  in  the  box  represents  the  mean. 

The  distance  from  the  lower  to  the  upper  quartile  is  called  the  interquartile 
distance,  and  it  will  be  represented  by  Q.  The  points  at  the  ends  of  the  two  lines 
(called  whiskers)  are  the  smallest  and  largest  points,  respectively,  within  1.5Q  of  the 
quantiles.  The  points  beyond  the  whiskers  are  outlying  values." 

Figures  1  and  2  present  the  boxplots  for  the  CD  hazards.  Figure  3  presents 
boxplots  for  the  neoplasm.  The  boxplots  are  grouped  by  level  of  TCE  which  is 
indicated  at  the  bottom  of  the  figure.  The  left  boxplot  in  each  grouping  is  for  the  3 
month  hazard.  The  right  boxplot  in  each  group  is  for  the  6  month  hazard. 
Comparison  of  the  boxplots  for  the  3  and  6  month  hazards  in  Figures  1  and  2 
suggests  that  the  3  and  6  month  hazards  are  roughly  the  same.  Comparing  Figures  1 
and  2  suggests  that  the  pretreatment  with  DEN  tends  to  increase  the  hazard. 
Comparison  of  the  3  and  6  month  boxplots  in  Figure  3  suggests  the  respective 
hazards  are  the  same  except  for  the  6  month  hazard  at  the  100%  TCE  level,  which 
appears  to  be  somewhat  higher  than  the  others  for  neoplasms. 

Figure  4  presents  the  histograms  of  the  CD  hazard  bootstrap  samples.  Once 
again  the  major  effect  seen  is  the  increase  in  hazard  for  the  fish  pretreated  with  DEN. 

Figure  5  presents  the  histograms  of  the  neoplasm  hazard  bootstrap  samples. 
Once  again  the  only  histogram  that  appears  different  is  the  histogram  for  the  6 
month  hazard  at  100%  TCE. 

Conclusions 

The  general  conclusion  from  the  above  analysis  is  that  there  is  only  a  weak  effect 
from  TCE  treatment  change,  regardless  of  whether  DEN  is  used.  The  effect  of  DEN  is 
noticeable:  the  second  3  months'  hazard  is  always  somewhat  larger  when  DEN  is 
used  than  is  the  case  with  no  DEN.  This  is  anticipated,  but  the  quantitative  degree  of 
enhancement  may  be  of  interest. 
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TABLE  3 

NONPARAMETRIC  HAZARD  FOR  CYSTIC  DEGENERATION 
(BOOTSTRAP  STANDARD  ERROR) 


Estimated  Hpard  _ 

(Standard  Error)  l^pq/2 5  J 

Sacrifice  Time 

Group 

DEN 

No  DEN 

%  TCE 

3  Months 

6  Months 

1 

X 

0 

0.24 

(0.09)  10.091 

0.04 

(0.11) 

3 

X 

25 

0.16 

(0.07)  [0.071 

0.27 

(0.16) 

5 

X 

50 

0.08 

(0.05)  [0.051 

0.22 

(0.14) 

7 

X 

100 

0.12 

(0.07)  [0.061 

0.11 

(0.11) 

2 

X 

0 

0.44 

(0.10)  [0.101 

0.11 

(0.20) 

4 

H 

25 

0.16 

(0.07)  [0  071 

0.54 

(0.17) 

6 

X 

50 

0.24 

(0.09)  [0.09] 

0.23 

(0.18) 

8 

X 

100 

0.28 

(0.09)  [0.091 

0.13 

(0.19) 

TABLE  4 

NONPARAMETRIC  HAZARD  FOR  NEOPLASMS 
(BOOTSTRAP  STANDARD  ERROR) 


Estimated  Hazard  1 
(Standard  Error)  \_^pqJ25  J 

Sacrifice  Time 

Group 

%  TCE 

3  Months 

6  Months 

2 

0 

0 

(0)101 

0.08 

(0.07) 

4 

25 

0.08 

(0.05)  [0.051 

0.16 

(0.12) 

6 

50 

0 

(0)  [0] 

0.08 

(0.07) 

8 

100 

0.08 

(0.06)  [0.051 

0.46 

(0.16) 

PARAMETRIC  ANALYSIS 


In  the  present  context  a  parametric  analysis  of  data  means  that  a  particular 
mathematical  form  is  adopted  for  the  distribution  of  T,  the  time  to  symptom 
occurrence.  It  is  desirable  that  such  a  form  have  a  plausible  biological  origin,  i.e.  that 
it  can  be  derived  from  suitable  biological  considerations,  and  that  it  adequately 
represent  the  data.  The  Moolgavkar  et  al.  models  (1973),  (1979),  (1983)  seem  to  satisfy 
the  former  requirement,  but  involve  at  least  four  parameters,  which  is  too  many  to 
attempt  to  fit  using  data  from  the  present  design.  Instead,  the  simple  exponential 
distribution. 


G(f;  A)  =  1  -r** 

has  been  adopted  for  illustration.  Note  that  the  single  parameter.  A,  is  actually 
interpretable  as  the  inverse  of  the  mean  of  T  (time  to  symptom  occurrence)  in  the 
population.  If  this  model  agrees  reasonably  well  with  the  data  then  1 /estimated 
A=l/Ais  easily  understood  and  interpreted.  The  exponential  model  also  implies 
that  the  theoretical  first  and  second  3  month  hazards  are  the  same.  Notice  that  since 
no  actual  times  to  symptom  appearance  are  ever  observed  such  a  quantity  is  not 
available  from  non-parametric  methodology.  The  parameter  A  (actually  it  is  best  to 
estimate  y=  log  A)  must  be  estimated  from  the  counts  at  three  months  and  six 
months.  The  method  used  here  is  that  of  maximum  likelihood ;  details  are  provided 
in  an  appendix. 

Tables  5  and  6  exhibit  the  results  of  the  analysis. 

These  results  seem  surprising,  since  mean  time  to  exhibit  the  CD  symptom 
appears  to  increase  with  TCE  dosage;  as  anticipated  the  effect  of  DEN  is  to  reduce  the 
time  to  symptom  appearance;  these  results  are  in  rough  qualitative  agreement  with 
the  non-parametric  results.  See  also  Figures  6-7,  which  indicate  the  uncertainty 
associated  with  the  above  numerical  values.  These  results  were  obtained  by 
bootstrapping. 

Figure  6  displays  boxplots  of  the  values  of  -y,  the  log  mean  time  to  CD  for  the 
bootstrap  samples.  Once  again  the  boxplots  are  grouped  in  pairs  by  level  of  TCE 
which  is  indicated  on  the  bottom  of  the  figure.  The  leftmost,  (respectively 
rightmost),  boxplot  in  a  group  is  for  the  fish  not  pretreated  with  DEN,  (respectively 
pretreated  with  DEN).  Once  again  the  major  effect  is  a  decrease  in  mean  time  to 
occurrence  of  CD  with  pretreatment  with  DEN.  The  variability  of  the  estimate 
makes  other  conclusions  suspect.  Figure  7  displays  the  histograms  of  the  bootstrap 
estimate  values  of  the  -y.  the  log  mean  time  to  CD. 

Figures  8-9  present  boxplots  comparing  the  bootstrap  estimates  of  the  probability 
that  CD  occurs  before  3  months  obtained  from  the  parametric  exponential  model  and 
the  nonparametric  analysis.  The  estimate  for  the  probability  using  the  exponential 
model  is 
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TABLE  5 

MAXIMUM  LIKELIHOOD  ESTIMATES  OF  MEAN  TIME  TO  EXHIBIT  CD 
(BOOTSTRAP  STANDARD  ERROR) 


Group 

DEN 

no  DEN 

%TCE 

Log  Mean  Time  to  CD 

CD  occurs  i 

[before  3  months  J 

(std  error) 

(std  error) 

1 

X 

0 

2.66 

(0.32) 

3 

X 

25 

2.68 

(0.34) 

0.19 

(0.05) 

5 

X 

50 

3.17 

(0.45) 

0.12 

(0.04) 

7 

X 

100 

3.19 

(0.60) 

0.12 

(0.04) 

2 

X 

0 

1.85 

(0.26) 

0.38 

(0.07) 

4 

n 

25 

23 

(0.31) 

0.26 

(0.07) 

6 

X 

50 

2.40 

(0.34) 

0.24 

(0.07) 

8 

X 

100 

2.33 

(0.33) 

0.25 

(0.07) 

TABLE  6 

MAXIMUM  LIKELIHOOD  ESTIMATES  OF  LOG  MEAN  TIME  TO  EXHIBIT 

NEOPLASMS 

(BOOTSTRAP  STANDARD  ERROR) 


EXPONENTIAL  MODEL 
PRETREATMENT  WITH  DEN 


Group 

%  TCE 

Log  Mean  Time  to 
Neoplasms 

1  Neoplasms  occurs 
[before  3  months 

2 

0 

4.97 

0.02 

(4.47) 

(0.02) 

4 

25 

3.34 

0.10 

(0.99) 

(0.04) 

6 

50 

4.97 

0.02 

(4.60) 

(0.02) 

8 

100 

2.88 

(0.39) 
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Pe  =1  -  expj-e^3j. 

The  estimate  of  the  probability  using  the  nonparametric  hazard  is  the  average  of  the 
first  and  second  3  month  hazard.  The  boxplots  are  grouped  by  level  of  TCE.  The  left 
(respectively  right)  one  in  each  group  are  the  bootstrap  estimates  for  the  parametric 
exponential  model  (respectively  the  nonparametric  hazard).  The  figures  suggest  that 
the  two  procedures  yield  roughly  the  same  estimate. 

Figure  10  presents  similar  boxplots  for  the  bootstrap  estimates  of  the  probability 
the  neoplasms  occur  before  3  months.  Note  that  the  exponential  model  estimates 
suggest  that  there  is  no  effect  at  100%  TCE.  Note  that  only  the  chronic  data  is  being 
examined. 

Figure  11  present  histograms  for  a  simulation  experiment  to  illustrate  the  effect 
of  using  more  fish  in  the  experiments.  Our  experiment  is  extreme  in  that  200  fish 
are  used  in  each  group;  100  are  sacrificed  at  3  months  and  100  sacrificed  at  6  months. 
The  nonparametric  estimates  of  G(3;0)  and  G(6;0)  for  each  group  of  the  CD  data  are 
used  as  the  true  probabilities  of  CD  occurring  at  3  and  6  months  respectively.  For 
each  simulation  replication  2  random  numbers  are  drawn;  one  from  a  binomial 
distribution  with  100  trials  and  probability  the  estimate  of  G  (3;0)  and  the  other  from 
a  binomial  distribution  with  100  trials  and  probability  the  estimate  of  G  (6;0).  For 
each  group  500  simulation  replications  are  done  and  the  two  3  month  hazards  are 
computed  for  each  replication  as  before.  A  comparison  of  the  histograms  in  Figures 
4  and  11  shows  the  amount  of  decrease  in  the  variability  of  the  estimates  that  can  be 
achieved  by  increasing  the  number  of  fish  used  in  the  experiment. 

BIOLOGICALLY-BASED  MODEL  DESCRIPTION 

It  is  widely  believed  that  pre-cancerous  conditions  in  an  organ  (the  liver)  occur 
as  a  result  of  cell  clonal  expansion,  followed  by  a  promotion  (to  tumor)  event. 
Specific  models  for  this  have  been  proposed  and  developed  by  Moolgavkar  and  co¬ 
workers.  More  recent  work  is  by  C.  J.  Portier  and  co-workers.  References  appear 
later. 

The  basic  mechanism  is  treated  as  random  or  probabilistic:  an  initiating  event, 
e.g.,  caused  by  contact  with  toxin,  affects  a  cell  within  an  organ  in  accordance  with  a 
simple  Poisson  process  with  rate  parameter  X.  That  is,  the  chance  of  an  uninitiated 
cell  being  initiated  in  time  interval  (f,  t+h )  is  approximately  Xh.  If  a  cell  is  initiated 
during  exposure  time  it  clones  itself  into  other  cells  at  rate  /3;  the  original  cells  and 
its  clones  die  randomly  at  rate  8.  All  cells  in  the  organ  perform  thus  independently, 
according  to  the  model.  Depending  upon  the  values  of  /?  and  8  (birth  and  death  rates 
respectively)  a  colony  of  initiated  cells  (pre-cancerous,  presumably)  either  tends  to 
grow  exponentially,  or  to  die  off  to  zero  (also  exponentially  fast).  The  fates  of 
colonies  characterized  by  the  same  values  of  birth  rate  and  death  rate  may  actually  be 
entirely  different,  as  befits  experience  with  variability  characteristic  of  the  real 
biological  world.  This  behavior  is  roughly  analogous  to  that  of  the  flipping  of  the 
same  coin:  on  one  occasion  10  flips  may  well  result  in  an  excess  of  5  Heads  (7  Heads 
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and  3  Tails),  analogous  to  more  births  (Heads)  than  deaths  (Tails);  on  another 
sequence  of  10  flips  with  the  same  coin  the  result  may  be  exactly  reversed  (7  Tails,  3 
Heads).  Processes  analogous  to  coin  flipping  or  dice  rolling  can  describe  much,  but 
possibly  not  all  interesting  biological  variability  pertinent  to  risk  analysis.  Other 
options  are  suggested  later. 

The  values  of  p  and  8  describe  clone  colony  properties  in  a  precise  probabilistic 
manner  if  the  model  is  correct.  It  is  certainly  only  approximate,  but  may  still  provide 
a  useful  tool  for  quantifying  risk  of  tumor  formation.  The  second  step  in  the 
malignant  cell  development  process  is  postulated  to  be  promotion.  A  model  for  this 
is  that  at  rate  p,  i.e.  with  probability  ph  in  time  (f,  t+h),  a  promotion  event  occurs  that 
affects  one  of  the  clone  colony  members  in  proportion  to  the  current  size  of  the 
colony;  such  events  are  assumed  to  occur  in  accordance  with  a  Poisson  process  with 
rate  proportional  to  instantaneous  clone  population  size.  At  the  instant  that  the  first 
such  promotion  event  occurs,  the  clone  colony  (if  one  exists,  i.e.  has  been  initiated) 
will  be  said  to  have  developed  a  tumor,  at  least  in  informal  layman's  terms.  Note 
that  all  original  cells  in  an  organ  are  assumed  to  be  independently  exposed  to 
initiation  and,  thereafter,  to  promotion.  Therefore  all  organ  cells  and  subsequent 
clones,  if  any,  must  survive  from  initiation  to  the  end  of  the  observation  period 
without  being  promoted  in  order  for  the  organ  to  survive  throughout. 

The  probabilistic  mechanism  described  has  been  used  to  obtain  a  formula  for  the 
survival  probability  for  an  organ  for  any  observation  time  t.  See  Appendix  B  for  the 
formula  and  its  derivation.  Similar  formulas  have  been  derived  also  by  Moolgavkar 
and  others.  Our  formula  provides  the  basis  for  statistically  estimating  from 
pathology  data,  (combinations  of)  the  parameters:  X,  the  initiation  rate;  p,  the 
promotion  rate;  and  P  and  8,  the  clonal  birth  and  death  rates.  Such  estimates  can,  in 
turn,  be  used  to  estimate  the  probability  of  cell,  and  organ,  survival  for  any  time 
period.  Appendix  A  contains  a  discussion  of  maximum  likelihood  estimation  from 
data  so  as  to  specify  parameters  of  a  preliminary  model.  Further  work  is  required  to 
obtain  additional  statistical  models  and  procedures  to  analyze  other  experimental 
data. 

Extensions  to  the  Model:  Extra-variation  of  Parameters 

The  above  model,  and  the  consequences  thereof  in  the  form  of  a  survival 
probability  function,  are  appealing  since  they  have  a  plausible  biological  basis. 
Organ-to-organ  outcomes  (tumor  occurrence  or  not)  vary  randomly,  but  according  to 
precisely  the  same  mechanism  in  each  organ;  i.e.  the  same  values  of  X,  p,  P  and  8  are 
assumed  to  hold  for  each  organ.  Note  that  this  ignores  likely  variability  between 
organs  in  different  subjects  (e  g.,  fish).  Different,  but  superficially  identical,  biological 
entities,  be  they  fish,  rats,  or  humans,  can  be  expected  to  have  some  differences;  these 
can  be  said  to  be  the  result  of  genetic  diversity.  Specifically  these  differences  may 
cause  the  effective  parameters  X,  pp  and  8  to  differ  substantially  across  animals.  If 
the  above  are  estimated  from  data  without  recognizing  the  possibility  of  extra¬ 
variation,  biased  results  will  be  obtained.  See  Harris  (1990)  for  biological 
explanations  of  inter-organ  (subject)  variability. 
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There  are  several  possible  simple  and  preliminary  ways  of  dealing  with  the 
above  problem.  One  is  by  attempting  to  "explain"  parameter  variation  by 
representing  it  as  a  function  of  some  causal  variable,  such  as  the  age,  sex,  weight,  etc., 
of  the  host  subject.  The  technique  is  a  variation  of  ordinary  regression  analysis; 
methods  of  McCullagh  and  Nelder  (1983)  suggest  themselves.  A  description  of  a 
preliminary  computational  procedure  to  estimate  model  parameters  is  described  in 
Appendix  A.  This  procedure  is  used  to  estimate  model  parameters  for  a  particular 
data  set.  A  second  approach  is  to  assume  that  the  variability  between  individual  host 
organs  can  be  represented  by  treating  some  or  all  of  the  parameters  as  random 
variables  with  their  own  distributions.  A  typical  survival  function  is  then  obtained 
by  mixing :  the  parametric  survival  function  of  Appendix  B  is  "simply"  randomized 
according  to  the  (joint)  distribution  of  the  parameters.  In  principle  it  is  desirable  to 
recognize  both  sources  of  variability  between  individuals,  adjusting  for  known 
sources  of  variation  by  a  regression  technique  where  possible,  but  recognizing  the 
"unexplainable"  variation  by  use  of  a  mixing  technique.  The  latter  has  been  carried 
out  to  a  limited  extent,  see  Gaver  and  Jacobs  (1992). 

CONCLUSIONS  AND  SUMMARY 

This  report  covers  an  initial  short  piece  of  research  conducted  under  the 
sponsorship  of  the  Army  Biomedical  Research  and  Development  Laboratory.  Its 
main  contribution  is  to  propose  and  illustrate  quantitative  assessments  of  treatment 
(here  groundwater  concentration)  effects  upon  medaka.  Those  quantifications 
include  the  estimation  of  statistical  sampling  errors  by  the  re-sampling  or 
bootstrapping  technique. 

The  somewhat  inconclusive  dose-response  relationships  revealed  seem  to 
imply  the  need  for  more  sensitive  experiments.  Possibly  such  sensitivity  can  be 
achieved  by  working  with  more  genetically  homogeneous  animals  (medaka). 
Possibly,  larger  numbers  of  animal  subjects  will  be  helpful  as  well.  Control  and 
measurement  of  experimental  conditions  (e.g.,  tank  temperature)  and  adjustments 
for  their  variations  can  play  a  useful  part  in  the  investigation. 

It  is  hoped  that  the  mathematical  and  statistical  approaches  illustrated  here  will 
help  to  promote  an  interest  in  the  further  use  of  such  ideas  among  biologists  and 
toxicologists. 
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APPENDIX  A 


MODEL  FITTING  METHODS  FOR  QUANTIFYING  BIO  ASSAYS 


PRELIMINARY  STATISTICAL  MODELS  AND  METHODS  FOR  ANALYZING 
BIOASSAY  DATA 

Suppose  N  organisms  (for  example  fish)  are  used  in  an  experiment.  Groups  of 
these  organisms  may  be  exposed  to  different  treatments.  Let  T(  be  the  random  time 
until  organism  i  develops  a  particular  symptom,  e.g.,  cystic  degeneration.  Let 
Xj  =  Xq,  ...,  X(-p)  be  covariates  which  (possibly)  influence  T,  ;the  X,-  could  be 
levels  of  substances  having  possible  toxic  effects  to  which  the  organisms  are  exposed. 
Let  G(f;xi)  =  P  |T(  <  f  | X,  =  Xj  J.  We  will  assume  that  the  organisms  develop  symptoms 

independently  of  each  other.  In  this  initial  model,  the  symptom  is  either  present  or 
not. 


Suppose  that  nk  organisms  are  sacrificed  at  time  tk  with  fj  <  f2  <  ...  <  tK  .  We 
will  label  the  organisms  so  that  organisms  1  through  n ^  are  sacrificed  at  time  f1 ; 
organisms  Hj+1, ...,  rij+w  2  are  sacrificed  at  time  f2;etc  Let  s,  =  1  if  organism  i 
exhibits  the  symptom  when  it  is  examined.  Under  the  assumption  of  independence, 
the  likelihood  function  is 

K  nk  s"jt-l+< 

1  =  IT  IT  ^ G(tk '  *nk- 1 +« ) 

Jt=l  i=l 

where  n0  =  0  and  G(t;  x)  =  1-G(f;  x).  The  likelihood  functions  form  the  basis  for 
estimation  of  parameters  in  the  distributions  that  model  survival  times,  i.e.  G. 

Example  (Simple  Binomial  Model).  If  there  are  no  covariates,  then  (A.l)  becomes 

L=ft  ("f^k),kC(tk)nt~,k  <A.2> 

where  fk  is  the  number  of  the  nk  organisms  exhibiting  the  symptom. 

A  procedure  to  estimate  the  parameters  of  the  distribution  G  for  the  simple 
binomial  model  is  as  follows. 


(A.l) 
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MAXIMUM  LIKELIHOOD  ESTIMATION  IN  THE  SIMPLE  BINOMIAL  MODEL 


(a)  Likelihood  and  Parameter  Estimation  Formulas 


Assume  the  distribution  of  the  time  to  appearance  of  a  symptom,  G,  is  a 
function  of  the  parameters  6j ,  j  -  1,  J.  In  this  section  we  discuss  maximum 
likelihood  estimation  of  Qj  for  the  simple  binomial  model.  Presumably  the  nk 
subjects  examined  at  time  tk,  k  =  1, 2, K  have  all  been  subjected  to  a  common 
dosage  of  a  potential  toxin.  The  purpose  of  the  present  analysis  is  to  predict  survival 
probabilities  as  they  depend  on  such  dosage.  The  log-likelihood  function  for  the 
simple  binomial  model  is 


K  f  . 

-  5>[  ?  +/tlnG((t;e)  +  (nt  -/*)taC( tt;9) 

*=1 


where  0  =  (0^, ..., 0j).  Differentiating,  we  obtain 


'■I 


dei  j£“G(fJt'e)a0/ 


G(i»;8) 


(A  .3) 


feh-Gfa-9)l-(”t-fe)Gfa;e) 

G(it;6)G((t;e) 


30, 


■GM) 


I 

/c=H 


fk-nkG(tk;9) 

3 

G(tk;0)C(tk;6) 

30y 

(A.4) 


Since  E[fk]  =  nkG{tk;9) 


E 


a2 

dOjde 


m 


^G(t,;e)^-G(t,;e) 

G{tk-,Q)G{tk;9) 


(A.5) 


Thus  a  Newton  procedure  for  finding  the  maximum  likelihood  estimates  of 
{0^;  j  =  1, ...,  /}  would  iteratively  solve  the  system  of  linear  equations 


O  =  -—/(0°)+  Y  E 
/  m=l 


30;30m 


(»»-«£) 


(A.6) 
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where  6°  =[d®,...,d(j j.  Such  iterative  procedures  can  be  programmed  for  a  digital 

computer,  and  the  resulting  parameter  values  can  be  used  to  compute  predictions  for 
survival  probabilities,  or  risk,  as  the  latter  depend  upon  the  parameters  of  such 
models  as  described  in  Appendix  B. 
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APPENDIX  B.  TWO-STAGE  CLONAL-EXPANSION  MODEL 


In  this  appendix  we  present  a  birth-death  model  for  the  distribution  of  time 
until  a  normal  cell  becomes  promoted  to  a  tumor. 

We  first  develop  an  expression  for  the  distribution  of  random  time,  S,  until  an 
initiated  cell  or  one  of  its  descendants  becomes  malignant. 

Assume  that  there  is  one  initiated  cell  at  time  0.  Such  cells  divide  at  an 
exponential  rate  /J,  and  die  at  an  exponential  rate  5.  Any  initiated  cell  turns 
malignant  at  an  exponential  rate  ji;  i.e.  is  the  promotion  rate. 

(a)  Time  to  Promotion  of  an  Initiated  Cell 

Let  S  be  the  random  time  at  which  some  initiated  cell  or  its  descendent  turns 
malignant;  note  that  S  may  actually  be  infinite  if  the  population  of  initiated  cell  and 
its  descendents  dies  out.  Put 


z(t)  =  P{S  >  t). 

The  following  probability  argument  provides  an  equation  for  z(f):  the  event  that 
S  >  t+A  {A  >  0)  occurs  if  (i)  neither  birth  (cloning),  death,  or  promotion  occurs  in 
(0 ,  A)  and  promotion  does  not  occur  in  (A,  t+A);  the  probability  of  this  is 
[1  -  (P+&+p)A  +  o(d)]z(f);  or  (ii)  birth/cloning  occurs  in  (0 ,A)  and  no  promotion  occurs 
in  {A,  t+A);  the  probability  of  this  event  is  [pA  +  o(d)]z2(f),  where  the  square 
recognizes  that  at  time  A  there  are  now  two  independent  clonal  families  to  be 
considered;  or  (iii)  the  original  initiated  cell  dies  in  (0,  A),  the  probability  of  which  is 
SA  +  o(A).  Sum  these  three  terms  to  obtain  the  probability  that  S  >  t+A: 

z(t  +  A)  =  (l-</J  +  5  +  /i)d)z(f)  +  p Az2(t)  +  8A. 

Now  subtract  z(f)  from  each  side,  divide  by  A  and  let  A  — >  0.  The  result  is  the 
differential  equation 


Hence  z(f)  satisfies  a  Riccati  equation  with  initial  condition 

z(0)  =  1 


(B.l) 


(B.2) 


The  solution  to  (B.l)  with  initial  condition  (B.2)  is 
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2(0  = 


Pl(l-P2)-P2(l-Plk/i(','-'’2>' 


where  pl  2  are  the  solutions  to  the  quadratic  equation 


(B.3) 


2 

x  -  1  +  -  +  — 

l  P  P 


H-* 


(BA) 


P\,2  ~  2 


i  5  p 
1  +  — +  £ 

P  P 


8  p 
1  +  — +  — 

P  P 


v  A 
J  P 


(B.5) 


Since 
p2  ^  1  and 


8  2  5"|1/2  8 

(l+|+|)  -4^  <  (l  +  ^+  jjj),  both  pT  and  p2  are  positive.  Further 


P\~P2  = 


fl  5  A5 

{  P  P)  P 


>0. 


Hence, 


Iim2(f)  =  p2.  (B.6) 

t—*  oo 

If  the  death  rate  5  =  0,  then  p,  =  0  and  lim  PfT  >  f}  =  0;  if  5  =  0,  then  there  is  no  death 

*  f  — »  O o 

of  initiated  cells  and  thus  an  initiated  cell  will  transition  to  a  malignant  cell  in  a 
finite  time  with  probability  1.  If  5>  0,  then  the  initiating  cells  can  die,  thus 
preventing  a  transition  to  malignancy  and  hence  lim  P(T >t)  =  p2> 0. 

(b)  Model  for  the  Time  until  a  Normal  Cell  becomes  Malignant  (is  Promoted  to 
Tumor) 

Assume  that  each  normal  cell  is  initiated  at  an  exponential  rate  Aq.  Let  N  be  the 
total  number  of  normal  cells  in  an  organ.  Let  T  denote  the  first  time  a  normal  cell 
transitions  to  a  malignant  cell. 
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P{T>f}  = 


i 

+  J*  Age~*Osz(f  -  s)ds 


N 


(B.7) 


where  z  is  a  given  in  (B.3).  Assume  Ag  is  small  and  put  A  =  AqN,  a  constant.  Then 

t 


P{T  >  f}  =  exp 


Nln 


0 


(B.8) 


=  exp 


t 


—At  +  A 


J  z(s)ds 


0 


(B.9) 


=  exp 


A(Pl-l)t-Ailn 


l_p2+(p1_1)^(Pi-P2)t 


P\~P2 


(B.10) 
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PRETREATED  WITH  DEN;CHRONIC  DEGENERATION 
NON  PARAMETRIC  HAZARD 


21 


Figure  2 
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